% This code loads directly the dbf files generated by the python files.

clear
% Base name of the grid:
BaseName = 'G1kmMiR';
% Year of reference for the IBGE geocodigo:
anorefer = 2007;
% Choose simplification of climate data:
SimpleOption = 2; % = 1 (Quarterly)
                  % = 2 (Harvest/Growth seasons)
                  % Definition of Harvest/Growth seasons in
                  % SimplifyGridClimateGrowthHarvest.m
regions = {'SP','CS'};

% Years of the sample:
Years{1}    = 2003:2013; % Years for SP
Years{2}    = 2004:2013; % Years for CS

for reg = [1,2]
    for y = Years{reg}
        % DBF file name.
        FileName                  = [BaseName, char(regions(reg)), num2str(y)];
        % Read DBF file.
        [DataTemp,FieldNamesTemp] = dbfread(FileName);
        % Convert variables of interest to table format.
        DataTemp                  = cell2table(DataTemp(:,[2,end]),'VariableNames', {'FID_Grid',['S',num2str(y)]});
        % Convert to dataset format.
        DataTemp                  = table2dataset(DataTemp);
        if y == min(Years{reg}) % for first year
            DS{reg} = DataTemp;
        else
            DS{reg} = join(DS{reg},DataTemp);
        end
        
        % Create identifiers for SP and CS
        DS{reg}.CanasatDataset = reg*ones(size(DS{reg}.FID_Grid));
        
        disp([' Sugarcane state loaded for ',char(regions(reg)),' year ', num2str(y)]);
    end
end
            
DS{2}.S2003 = NaN(size(DS{2}.S2004));

DS = [DS{1};DS{2}];


clear DataTemp FieldNamesTemp
% Import Climate Micro data and simplify it:
if SimpleOption == 1
    SimplifyGridClimateQuarterly;
else
    SimplifyGridClimateGrowthHarvest;
end
disp(' ')

% Import Grid to Municipality table:
disp(' Importing Grids to Municipalities ')
GridMunic = dbfread([BaseName,'Munic']);
GridMunic = cell2table(GridMunic,'VariableNames', {'Join_Count','FID_Grid','JOIN_FID','Join_Count2','TARGET_FID','anoderefer','munic'});
GridMunic = table2dataset(GridMunic);
GridMunic = GridMunic(GridMunic.anoderefer == anorefer,:);
disp(' ')

% Import Grid to Sugarcane Road Distance table:
disp(' Importing Grids to Sugarcane buffers ')
[GridSugRoadDist,NamesSugRoadDist] = dbfread([BaseName,'SugRoadDist']);
NamesSugRoadDist{2} = 'FID_Grid'; % Give std name for Grid ID
GridSugRoadDist = cell2table(GridSugRoadDist,'VariableNames', NamesSugRoadDist);
GridSugRoadDist = table2dataset(GridSugRoadDist);
GridSugRoadDist(:,1) = []; % Eliminate first column, that adds no information. 

% Import Grid to Elevation:
disp(' Importing Grids to Elevation ')
GridElev = dbfread([BaseName,'WGSElevation']);
GridElev = cell2table(GridElev,'VariableNames', {'Join_Count', 'TARGET_ID', 'FID_Grid', 'Altitude', 'Slope'});
GridElev = table2dataset(GridElev);
disp(' ')

% Import Grid to GAEZ yield information:
disp(' Importing Grids to GAEZ yield ')
[GridGAEZ] = dbfread(['G1kmG5kmWGSGAEZ']);
GridGAEZ = cell2table(GridGAEZ,'VariableNames', {'G5km_FID', 'FID_Grid', 'gaez_aclim_mze', 'gaez_aclim_soy', 'gaez_aclim_suc', ...
    'gaez_aecol_mze', 'gaez_aecol_soy', 'gaez_aecol_suc', 'gaez_act_mze', 'gaez_act_soy', 'gaez_act_suc'});
GridGAEZ = table2dataset(GridGAEZ);
disp(' ')

% % Import Grid to Slope:
% disp(' Importing Grids to Slope ')
% GridSlope = dbfread([BaseName,'Slope']);
% GridSlope = cell2table(GridSlope,'VariableNames', {'FID_Grid', 'FID_Filter', 'FID_Estado', 'FID_Buffer', 'Slope'});
% GridSlope = table2dataset(GridSlope);
% disp(' ')

% Import Grid to Soil:
% disp(' Importing Grids to Soil ')
% GridSoil = dbfread([BaseName,'WGSSoil']);
% GridSoil = cell2table(GridSoil,'VariableNames', {'Join_Count', 'TARGET_ID', 'FID_Grid', 'CEC', 'CLYPPT', 'ORCDRC', 'PHIHOX', 'SLTPPT', 'SNDPPT'});
% GridSoil = table2dataset(GridSoil);
% GridSoil = GridSoil(:,{'FID_Grid', 'CEC', 'CLYPPT', 'ORCDRC', 'PHIHOX', 'SLTPPT', 'SNDPPT'});
% disp(' ')

% Import Grid to Nasa table:
% disp(' Importing Grids to Nasa ')
% [GridNasa,VarNames] = dbfread([BaseName,'Nasa']);
% GridNasa            = cell2table(GridNasa,'VariableNames', {'Join_Count', 'FID_Grid', 'nasacropland2000', 'nasapasture2000'});
% GridNasa            = table2dataset(GridNasa);
% GridNasa            = GridNasa(:,{'FID_Grid', 'nasacropland2000', 'nasapasture2000'});
% disp(' ')

% Import Grid to Urban table:
disp(' Importing Grids to Lakes ')
[GridUrban,VarNames] = dbfread([BaseName,'Urban']);
GridUrban            = cell2table(GridUrban,'VariableNames', {'Join_Count', 'FID_Grid', 'JOIN_FID', 'Urban', 'TARGET_FID'});
GridUrban            = table2dataset(GridUrban);
disp(' ')

% Import Grid to Lakes table:
disp(' Importing Grids to Lakes ')
[GridLakes,VarNames] = dbfread([BaseName,'Lake']);
GridLakes            = cell2table(GridLakes,'VariableNames', {'Join_Count', 'FID_Grid', 'JOIN_FID', 'Lake', 'TARGET_FID'});
GridLakes            = table2dataset(GridLakes);
disp(' ')

%Import Grid to Carbon Stock table:
% disp(' Importing Grids to Carbon Stock ')
% [GridCarbonStock,VarNames]  = dbfread([BaseName,'WGSCarbonStock']);
% GridCarbonStock             = cell2table(GridCarbonStock,'VariableNames', {'Join_Count', 'TARGET_FID', 'FID_Grid', 'CarbonStock'});
% GridCarbonStock             = table2dataset(GridCarbonStock);

%Import Grid to Last of the Wild table:
% disp(' Importing Grids to Last of the Wild ')
% [GridLastOfWild,VarNames]  = dbfread([BaseName,'LastOfWild']);
% GridLastOfWild             = cell2table(GridLastOfWild,'VariableNames', {'Join_Count', 'FID_Grid', 'JOIN_FID', 'ECO_NUM', 'G200_NUM', 'Join_Cou_1','TARGET_F_1'});
% GridLastOfWild             = table2dataset(GridLastOfWild);
% GridLastOfWild             = GridLastOfWild(:,{'FID_Grid','ECO_NUM','G200_NUM'});

% Import Grid to 40km Buffer table:
disp(' Importing Grids to 40km Buffer ')
[Grid40kmBuffer]  = dbfread([BaseName,'Buffer40km']);
Grid40kmBuffer             = cell2table(Grid40kmBuffer,'VariableNames', {'Join_Count', 'FID_Grid', 'TARGET_FID' , 'Buffer40km', 'fid'});
Grid40kmBuffer             = table2dataset(Grid40kmBuffer);
Grid40kmBuffer             = Grid40kmBuffer(:,{'FID_Grid','Buffer40km'});

% Import Grid to Cost Distance Def1 table:
% disp(' Importing Grids to Cost Distance Def1 ')
% [GridCostDistanceDef1,NameCostDistanceDef1] = dbfread([BaseName,'WGSCostDistanceDef1']);
% GridCostDistanceDef1           = cell2table(GridCostDistanceDef1,'VariableNames', {'Join_Count', 'TARGET_FID', 'FID_Grid', 'CostDistanceDef1'});
% GridCostDistanceDef1           = table2dataset(GridCostDistanceDef1);
% GridCostDistanceDef1 = GridCostDistanceDef1(:,{'FID_Grid' 'CostDistanceDef1'});

% Import Grid to Road Cost Distance:
disp(' Importing Grids to Cost Distance Def1 ')
[GridRoadCostDistance,NameCostDistanceDef1] = dbfread([BaseName,'WGSRoadCostDistance']);
GridRoadCostDistance           = cell2table(GridRoadCostDistance,'VariableNames', {'Join_Count', 'TARGET_FID', 'FID_Grid', 'RoadCostDistance'});
GridRoadCostDistance           = table2dataset(GridRoadCostDistance);
GridRoadCostDistance = GridRoadCostDistance(:,{'FID_Grid' 'RoadCostDistance'});

% disp(' Importing Grids to road distance and transport cost to ports ')
% GridRoadCostDistanceAll = readtable('G1kmMiR_TransCost2Port.csv');
% GridRoadCostDistanceAll.Properties.VariableNames{'TARGET_FID'}   = 'FID_Grid';
% GridRoadCostDistanceAll.Properties.VariableNames{'RoadDistance'} = 'RoadDist2Port';
% GridRoadCostDistanceAll = GridRoadCostDistanceAll(:,2:end);
% GridRoadCostDistanceAll = table2dataset(GridRoadCostDistanceAll);
% [GridRoadCostDistanceAll,NameCostDistanceDef1] = dbfread([BaseName,'_TransCost2Port']);
% GridRoadCostDistanceAll           = cell2table(GridRoadCostDistanceAll,'VariableNames', {'FID_Grid','RoadDist2Port', 'Corn2Port', 'Soy2Port', 'Sugar2Port'});

disp(' Importing Grids to protected areas ')
GridProtected = readtable('G1kmMiRProtected.csv');
GridProtected.Properties.VariableNames{'TARGET_FID'}   = 'FID_Grid';
GridProtected = table2dataset(GridProtected);

disp(' Importing Grids to MapBiomas ')
MapBiomas = readtable('G1kmMiRMapbiomas.csv');
MapBiomas.Properties.VariableNames{'TARGET_FID'}   = 'FID_Grid';
MapBiomas = table2dataset(MapBiomas);


% Import Grid to Coordinates:
disp(' Importing Grids to Coordinates ')
GridXY = dbfread([BaseName,'XY']);
GridXY = cell2table(GridXY,'VariableNames', {'Join_Count', 'FID_Grid', 'CoordX', 'CoordY'});
GridXY = table2dataset(GridXY);
GridXY = GridXY(:,{'FID_Grid','CoordX','CoordY'});
disp(' ')

disp(' ')

disp(' Creating GridInfo with additional information on grids: ')
disp(' Urban Areas, Lakes, Municipalities, Elevation, Slope, Carbon Stock, Transportation cost, 40 km Buffer, ')
disp(' and Soil')
GridInfo = join(GridMunic(:,{'FID_Grid','munic'}),GridUrban(:,{'FID_Grid','Urban'}));
GridInfo = join(GridInfo,GridSugRoadDist);
GridInfo = join(GridInfo,GridLakes(:,{'FID_Grid','Lake'}));
GridInfo = join(GridInfo,GridCostDistanceDef1);
GridInfo = join(GridInfo,GridRoadCostDistance);
GridInfo = join(GridInfo,GridRoadCostDistanceAll);
GridInfo = join(GridInfo,GridProtected);
GridInfo = join(GridInfo,MapBiomas);
GridInfo = join(GridInfo,GridElev(:,{'FID_Grid','Altitude','Slope'}));
GridInfo = join(GridInfo,GridCarbonStock(:,{'FID_Grid','CarbonStock'}));
GridInfo = join(GridInfo,Grid40kmBuffer);
GridInfo = join(GridInfo,GridNasa);
GridInfo = join(GridInfo,GridLastOfWild);
GridInfo = join(GridInfo,GridClimate);
GridInfo = join(GridInfo,GridGAEZ);
GridInfo = join(GridInfo,GridSoil);
GridInfo = join(GridInfo,GridXY);
disp(' ')

disp(' Join GridInfo to sugarcane state yearly dataset ')
try
    DS = join(DS,GridInfo,'FID_Grid');
catch
    disp(' Disparity between Field ID in GridMunic and DS ')
    disp(' I found this problem when using a very thin grid (500m).') 
    disp(' One grid point was exactly at the border of two municipalities ')
    disp(' Know ID this point and eliminate ')
    [C,ia] = setdiff(DS(:,'FID_Grid'),GridInfo(:,'FID_Grid'),'FID_Grid')
    
    ind2keep = (DS.FID_Grid>0);
    for in=ia'
        ind2keep = (ind2keep & (DS.FID_Grid~=DS.FID_Grid(in)));
    end
    DS = DS(ind2keep,:);
    
    DS = join(DS,GridInfo,'FID_Grid');
end



% Save Dataset Properties:

% Choose simplification of climate data:
if SimpleOption == 1; % = 1 (Quarterly)
    Simplification = 'quaterly';
else
    Simplification = 'harvest/growth seasons';% = 2 (Harvest/Growth seasons)
end

Description = ['CS, Grid: ' BaseName ', Years: ' num2str(min(Years{1})) ' to ' num2str(max(Years{2})) ', munic reference: IBGE geocodigo year ' num2str(anorefer) ', Climate data: ' Simplification];
DS.Properties.Description = Description;
DS = DS(:,[1 3 2 4:end]); % Rearrange columns so that cana state var are kept all together 
save(BaseName,'DS','-v7.3');